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Abstract 



2^ ■ The fluctuation exchange (FLEX) approximation is apphed to study the 

■ Holstein-Hubbard model. Due to the retarded nature of the phonon-mediated 

electron-electron interaction, neither fast Fourier transform (FFT) nor previ- 
^ ■ ously developed NRG methods for Hubbard-type purely electronic models are 



applicable, while brute force solutions are limited by the demands on com- 



o 

O ■ putational time and storage which increase rapidly at low temperature T. 



Here, we describe a new numerical renormalization group (NRG) technique to 



■ solve the FLEX equations efficiently. Several orders of magnitude of CPU 

time and storage can be saved at low T (~ SO-fT). To test our approach, 
we compare our NRG results to brute force calculations on small lattices at 
elevated temperatures. Both s-wave and d-wave superconducting phase di- 
agrams are then obtained by applying the NRG approach at low T. The 
isotope effect for s-wave pairing is BCS-like in a realistic phonon frequency 
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range, but vanishes at unphysically large phonon frequency (~ band width). 
For d-wave pairing, the isotope exponent is negative and small compared to 
the typical observed values in non-optimally doped cuprates. 
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I. INTRODUCTION 



A large body of experimental evidence suggests that several liigh-Tc cuprate super- 
conductors exhibit a pairing state of dx2-y2 symmetryllli. In combination with the large 
superconducting transition temperature of these materials, this suggests a superconduct- 
ing pairing mechanism of predominantly electronic origin. In particular anti-ferromagnetic 
(AF) spin fluctuation (SF) exchange has been proposed as a possible electronic candidate 
mechanismim which would tend to give rise to (1^2 _y2 pairing symmetry. 

However, except near certain "optimal" doping concentrations, many cuprates ex- 
hibit a quite noticeable doping dependent isotope effectlH and other pronounced, 
superconductivity-related lattice effect^!. This indicates that electron-phonon (EP) in- 
teractions could be important and should be included in the theory.0'lll 

In the past decade, conserving self-consistent field (SCF) methodslll and related dia- 
grammatic approaches, such as the fluctuation exchange (FLEX) approximationilMUHl have 
been used to study AF SF exchange within the framework of microscopic correlated electron 
models. Most of the FLEX-based SCF studies so far have been limited to Hubbard-type 
models with instantaneous, local Coulomb interactions in simple tight-binding models, and 
some extensions to long-range Coulomb interactionsil'ii. For these types of model systems, 
both numerical renormalization group (NRG)§ and fast Fourier transform (FFT)li meth- 
ods have been developed and successfully applied to solve the FLEX equations efficiently on 
large space-time lattices in the physically relevant low-temperature regime. The numerical 
solution of the FLEX equations is greatly simplified in instantaneous interaction models, due 
to the fact that the bare electron-electron interaction potential is frequency independent. 

When EP interactions are introduced into the model, already the bare electron-electron 
potential becomes explicitly frequency dependent, due the retarded nature of the phonon- 
mediated interaction. In that case, the solution of the FLEX equations requires the inversion 
of certain, large fermion frequency matrices which, in a brute force approach, would increase 
the demands on CPU time and memory consumption by several orders of magnitude, relative 
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to the purely electronic models studied so far. Also, neither FFT approachesllZl nor the 
original form of the NRG methodEl can be used here to reduce the numerical effort to 
manageable proportions. 

Previous SCF studies of strongly correlated electron models with EP couplingil have 
included renormalization of the Coulomb and phonon-mediated potential at the level of an 
RPA particle-hole bubble. In this much simpler approach, one neglects the electron-electron 
exchange scattering which arises in the full FLEX approximation, due to the Pauli exclu- 
sion principle. This simpler RPA-based approximationii then avoids the frequency matrix 
inversion problem, since the latter arises only from the exchange ladder diagrams. This may 
be a good approximation when the Coulomb repulsion is strongly screened and the phonon 
frequency is small. These are, essentially, the conditions under which Migdal's theorem is 
validii and the Eliashberg theoryiHH§ of phonon-mediated conventional superconductivity is 
applicable. In general, these conditions may not be satisfied in correlated electron systems 
and it may be necessary to include both Coulomb and EP contributions to the electron 
exchange vertex. 

In the present paper, we describe an extension of the NRG approach which will allow us 
to incorporate EP interactions into a Hubbard-type correlated electron model and handle 
the resulting frequency matrix inversions in the FLEX equations efficiently. An efficient 
algorithm to solve this matrix inversion problem in the present context is also an important 
first step towards studying the next level of SCF theory, such as, for example, the parquet 
theoryiiiil. Our present treatment is limited to the case of the so-called Holstein-Hubbard 
modelii'ii'0 where the phonon-mediated electron-electron potential is without momentum 
dependence. However, when combined with recently developed real-space basis representa- 
tion approaches,ii'00 our basic method to the frequency matrix inversion problem will also 
become applicable to bare potentials which are frequency and momentum dependent. 

Preliminary results obtained with our present frequency matrix NRG technique for the 
(i-wave instability of the Holstein-Hubbard model have been reported elsewhereil. The 
purpose of the present paper is to give a first detailed account of the technique itself, to 



present results for the s-wave instability and for the competition between d- and s-wave 
pairing in their respective parameter regimes. 

The paper is organized as follows: in Section II, we summarize the FLEX approximation 
for the Holstein-Hubbard model and define the notations used in the paper. In Section 
III, we describe our fermion frequency matrix NRG technique in detail. In Section IV, we 
present our results obtained with the NRG approach for the Holstein-Hubbard model. They 
include a comparison of results for one-particle correlation functions obtained with the NRG 
and, respectively, by the brute force approach; and some applications to the superconducting 
instabilities and the isotope effect of the Holstein-Hubbard model. We conclude with a brief 
summary in Section V. 

II. HOLSTEIN HUBBARD MODEL IN THE FLEX APPROXIMATION 

We start from the simplest microscopic model, including both an on-site Hubbard U 
Coulomb repulsion and a local Holstein-type EP coupling to an Einstein phonon branch. 
The Hamiltonian of this Holstein-Hubbard modeliiiiii can be written as: 



with a nearest neighbor hopping t, chemical potential /x, on-site Coulomb repulsion U, on- 
site EP coupling constant C, harmonic restoring force constant K, and ionic oscillator mass 
M. The cj^ (Cj^) is the electron creation (annihilation) operator at site i and spin a; rii^ is the 
number operator; and Ui is the local ionic displacement at lattice site i and Pi = —ihd/dui. 
The dispersionless bare phonon frequency is 




(1) 



VLp = (K/M) 



1/2 



(2) 



and the phonon-mediated on-site attraction is 



Up = C^/K . 



(3) 
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The bare interaction vertices entering into the FLEX treatment are shown in Fig.|l|, 
including the particle-hole density and magnetic vertices [ Figs, [^(b) and 0(c)] 

^nx,n2,n3,n4i'^^rn) = '^^m ,714 (*'^m)'5ni,n2+m^n3+m,n4 , (4) 

where 
and 

^ni,n2,n3,n4i'^^rn) = '^ni,n4 (^'^ni)'^™i,"2+rra'^n3+m,n4 ; (6) 

where 



V. 



nuaS^^m) = -bp(^^ni " ^^n^) + U] , (7) 



and the particle-particle singlet and triplet vertices [ Fig. |T](d)] 

^i,"2,"3,"4(^^'^) ~ '^ni,n4i'^^iri)Sni,-n2+m^-n3+m,n4 1 (8) 

where 
and 

^l,n2,n3,n4(^^»Ti) ~ '^ni,n4 (^'^"i)'^'^i,-"2+m'^-n3+m,n4 (10) 

where 

^ni,n4(^i^m) = ^[t;p(«t^„,i - iu;„J - t;p(ia;„i + - f^m)] ■ (11) 

Here, the electron-electron potential includes both the the Hubbard U and the phonon- 
mediated contribution [ Fig. |l](a)] 

vpiium) = -Upnl/inl + ul) . (12) 



The boson Matsubara frequencies are denoted by Vm = ^muT and the fermion Matsubara 
frequencies by u;„ = {2n + 1)'kT with integer m and n. 

In the FLEX approximation, the single-particle self-energy is then given byB 



T 



+VPP{k-k';iuJn)G*{k')} 



(13) 



where the effective interaction potentials are given b; 



V2{q;iujr. 



-Vp{iUrn) + 



V^^^{q- luj^) = - E { [^(1 + 5)"' - S\ Iq) vl„^^{%v^) + 
3[r(l + T)-1-T] ,(g) t;J,,„(2z/j} , 
Xp/,(g; zu;n") 



(15) 



Rn,n"{q) = V „{Wm) 



(16) 
(17) 



r D or M , 
for = <^ 

Here, we are using a momentum-energy vector notation where, for fermion lines, k = 
(k,iujn) and k' = {k.',iujn') and, for boson lines, q = (q, «//„). The Green's function 
is 



G{k) = [iUn-e^-J^ik)]-' 



and the tight binding band 



ek = — 2t(cos + cos k^,) — /i . 
The bare particle-hole and particle-particle fluctuation functions are then defined as: 



Xph{q;iuJn) 



Zj2Gik + q)Gik) 
k 

T 



N 



Y,G{k + q)G{-k) 



(18) 

(19) 

(20) 
(21) 



without summation over the fermion Matsubara frequency iuon- 

Eqs. (p!3|-pTD are solved iteratively. The iteration proceeds as follows: (1) choose the 



temperature and either a fixed electron concentration (n) or a fixed chemical potential /i; 
(2) guess an initial self-energy at this temperature and electron concentration, e.g., 
S = 0; (3) calculate the Green's function G{k) from Eq. (JT^); (4) calculate the bare 
particle-hole and particle-particle fluctuation functions by equations (^) and (|2ll) using 
the Green's function obtained in step (3); (5) evaluate the matrices R in all four channels 
and the effective potentials V2, V"^^, and V^^] (6) update the self-energy by Eq. (|l^); (7) 
using the updated self-energy S from step (6) as input, go back to step (3) to calculate an 
updated Green's function G. The iterative cycle, consisting of steps (3) through (7), is then 
repeated until a converged self-energy is obtained. If {n) is to be fixed to a given input 
value, then the chemical potential n must be aj dusted accordingly during the self-consistent 
calculation. 

Because of the retarded nature of Vp^iUm) the bare vertices in Fig.l(b-d) depend explic- 
itly on the internal frequency transfer. As a consequence, a frequency matrix inversion is 
necessary to evaluate the fluctuation potentials, V^^^ and V^^, in Eqs. (|TB|) and (P^. In a 
brute force approach, this matrix dimension can become as large as 500^ to 1000^ (the size 
of the entire fermion Matsubara frequency set) near the transition temperature. Further- 
more, for each iteration, the number of such matrix inversions to be carried out is about the 
same as the number of boson Matsubara frequencies times the size of the momentum grid. 
At the space-time lattice sizes required to study the physically interesting low-temperature 
regime, it is therefore not possible to carry out such a brute force calculation with currently 
available computing resources. Recent FFT and numerical NRG techniques, developed for 
the pure Hubbard FLEX equations are also not directly applicable. In the next section, we 
will describe a generalized "fermion frequency matrix" NRG method to handle the numerics 
efficiently. 

After a convergent self-energy is obtained at a fixed temperature and electron filling, 
the search for pairing instabilities requires the calculation of the maximal eigenvalue of the 
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pairing kernel, A(T), as a function of temperature T, from 

X{T)m = E V^,Uk.k';T)G{k')G{-k')(j){k') , (22) 

where the pairing potential V^air is^: 

T^pair(fc,fc';T) = ^v^(^u;^-^uJ^,)-^[D{l + D)-X^,ik-k')v^,^__,Xk-k') 

+ ^[M(l + MrXn'ik - k'Hl_„{k - k') . (23) 

The instability is reached when A(T) approaches unity, i.e., 

A(T) ^1 T ^ Tc. (24) 



III. NUMERICAL RENORMALIZATION GROUP APPROACH FOR THE 

HOLSTEIN HUBBARD MODEL 

A numerical NRG method has been successfully applied to the FLEX equations of the 
Hubbard mo del0, which has a frequency- and momentum-indepedent bare interaction, the 
on-site Coulomb U. In this case the matrix inversion in equations (|T5|) and (^) can be 
carried out analytically and the NRG operations are greatly simplifiedEl. Due to the fre- 
quency dependence of the phonon-mediated interaction Vp, we now have to construct a more 
general NRG operation in which the frequency dependence of the bare interaction is taken 
into account. The detailed procedure will now be described. 

The NRG evaluation of the self-energy follows closely the original NRG approach de- 
scribed in Ref. O. We will largely adopt the notation and terminology introduced therein. 



We are implementing a pure "frequency NRG" (in the sense of Ref. |T^). That is, the grid 
of momentum points k is chosen from the outset to be dense enough for the lowest temper- 
atures to be reached and remains constant throughout all NRG steps; only the Matsubara 
frequency grids {iun, ium) change from one NRG step to the next. 
The basic assumption underlying the NRG approach is that 
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(1) quantities such as the and G{k) are, to good approximation, independent of 
temperature at high frequencies, \iUn\'^T and 

(2) within that high-frequency regime, they are slowly varying with iuJn, on u;„-scales of 
order T and that 

(3) the contribution to ^{k) arising from scattering into the high frequency region 
[|ic<j„/|^T in Eq. ([13|) ], denoted by AT,{k) below, is to good approximation independent of 
temperature and slowly varying for all iun- 

For the case of the pure Hubbard model, these NRG assumptions have been verified in 
great detail, by explicit numerical calcuationsEl. A general justification of these assumptions 
can be given. It is based on the notion that the energy denominators in the Green's function 
G{k) become very large and essentially T-independent for \iuj \ ^T. Hence, all strongly T- 
dependent details are "washed out" in the high-frequency regimelli. 

As a consequence, only the low-frequency part of E, for \iujn \ ^T, arising from scattering 
into the low-frequency region \iun'\ ^ T in Eq. exhibits substantial T-dependence. 

As described in Ref. |TU|, one therefore divides the self-energy S(A;) in Eq. (|TB|) into two 
contributions, arising respectively from the scattering iojn^ioJn' into a "low" region L [i.e. 
iujni G L in equation (|TB|)] and into a "high" region [iuin' E H in equation (|r^)], that is, 

m = §E E S^ik,k')+ AE{k) . (25) 

k' iuj„i£L 



Here S^{k,k') denotes the summand in Eq. ([13D , and AS is the contribution from the 
iuJn' -summation over the "high" region H. 

The basic idea of the NRG approach is to reduce the numerical effort by evaluating 
AE(/c) at a higher temperature, on a correspondingly coarser Matsubara grid. This higher- 
T result is then interpolated onto the finer Matsubara grid relevant to the lower T. The 
interpolation onto the lower-T Matsubara grid needs to be performed only for Matsubara 
frequencies iun G L. Only the "low" contribution in Eq. ( pS]) needs to be evaluated by 
summing iun' over the L-portion of the finer, lower-T grid and this, again, needs to be done 
only for iUn^L. 
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Starting from a large initial temperature Tq and large initial Matsubara summation cut- 
off Qq, this basic NRG step is carried out repeatedly, through a sequence of decreasing 
temperatures 

To > Ti > ... > T, > ... (26) 

and decreasing Matsubara cut-offs 

no>ni> ... >n^> ... , (27) 

until the desired final temperature is reached. The initial maximal cutoff must be chosen 
large enough that the physical results of the calculation, T) and A(T), are independent 
of Qq, i.e., typically large compared to the bandwidth 8t. 

The subsequent renormalized fli (with i > 0) delineate the boundaries between the low 
and high regions, Lj and Hi, in the i-th NRG step where 

Li = {icu^^^ I ni > n integer} (28) 

and 

Hi = \J^Hj (29) 

i=i 

comprises the high-frequency region increments Aifj of the present (j = i) and all prior 
(j < i) NRG steps, given by 

^Hj = I ^-^i > \uj^J~^'^\ > n^; n integer} (30) 

The respective fermion Matsubara frequency grids are given by 

iu^i'^ = {2n + l)'KTj . (31) 

for integer n and j. In order to ensure that the summands in Eq. ( |55D below enter with the 
correct weights, the Qj and Tj are chosen to obey the grid matching conditions 

= 27: NjTj, Vlj = 2'KKjTj^i (32) 
11 



so that 



Tj/Tj_i = Kj/Nj, ^/^.i = Kj/Nj_i (33) 

with Nj and Kj integer and Kj > Nj > 1, a.s described in detail in Ref. |19[ The calculations 
presented below are based on a "factor-2" renormalization group with Qj/Qj_i = Tj/Tj_i = 
1/2, as illustrated by the frequency grids shown in Fig. ^ Note that the resulting NRG 
fermion Matsubara grid Li U Hi at temperature Tj is substantially "diluted" compared to 
the full, dense Matsubara grid 

Ai = I fio > kn^l; n integer} . (34) 

The latter constitutes the basic frequency summation domain in a "brute force" calculation 
at temperature Tj. 

The contribution from the Hi region, AS*^*^(A;), is "frozen in" during the self-consistency 
iteration at temperature Tj. It is calculated, prior to the self-consistency iteration, by 

AE»(A;) = g^E E S^\k,k') (35) 

j=0 k' iJ^)^^H,+, 

n' 

= AS(-^) + %i^ E St'\hk') (36) 

n' 

= AS(*-i) + , (37) 

where now k' = (k', ict;^^?) and S^\k,k') is the summand in Eq. (pTSD, computed with the 
Green's function G^^\k'), which, in turn, is obtained via the Dyson equation (|18D from the 
self-energy T^^^^k') on the grid Lj UHj at temperature Tj. Note here that the linear interpo- 
lation of the summand from higher- to lower-T frequency grids introduces the temperature 
re- weighting factor Tj/Ti into the summation in Eq. (|35D. 

The evaluation of the increment ^S*^*-* via Eq. ( pTf ) needs to be carried out only for iu- 
points in Both and S'-*^^-' are then added and interpolated onto the finer iu^^^-grid, 
inside Li, to obtain AS*^*-* on Li. This fixed AE^*) is then used in the (re-) calculation of S^*) 
on Li, via Eq. (^Sf), during the self-consistency iteration at temperature Tj. 
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Note that AS*^'^(k, iuj^"^) does not have to be calculated or stored for grid points outside 
of the low-frequency domain Lj, that is, for icol^'^ G Hi. The values of S(*)(k, itol^^) on Hi are 
already "frozen in" before the self-consistency iteration at temperatuure T,, since, by the 
above-stated NRG assumption of T-independece in Hi, 

E»(k,z4^')) = S(^-i)(k,z^(^')) = ... = S(^)(k,z4^')) (38) 

for i> j and ito^'' E Hi. Hence, the values of S*^-'^ on Hj C Hi become "frozen in" (i.e., 
permanently stored) after the j-th NRG step, and need not be re-calculated during sub- 
sequent NRG steps i> j. Only the values of E^*) on the low-frequency grid Li need to be 
re-calculated and iterated to self-consistency during the i-th NRG step. 

In order to evaluate the self-energy S, we therefore need the effective potentials V2, V^^ 
and V^^ at iuJn and iujn' in the L-region only. However, the fermion frequency summation 
over iUn" and the fermion frequency matrix inversions in equations (jl^) and (^6|) still run 
over the whole Matsubara frequency range up to the cutoff f^o- We therefore have to develop 
an efficient algorithm to overcome the fast growth of CPU time and memory requirements 
associated with these matrix operations at low temperatures. To this end, the NRG ap- 
proach, outlined above for the self-energy E(/c), must be extended to the fermion frequency 
matrices, i?„^„/(g), where R stands for D, M, S, or T, as defined in Eq. (P^). In the following, 
we will, for simplicity, omit the g-argument and use the notation 

R{iuJn,iuJn') = R{iuJn,iuJn'; O) = Rn,n'{q) ■ (39) 

As defined in Eqs. (|T7| - |21|) , these i?-matrix elements can be regarded as the values of 
an analytical function R{iuj,iuj'), defined for continuous iuo- and icu'-arguments. The basic 
assumption underlying our NRG approach for the i?-matrix is that ^R{iuj,iuj'), as defined 
in Eqs. (|l3)-(^ is 

(1) independent of temperature and 

(2) slowly varying on an iu-scale of order T, 

if either \iuj\ or \iuj'\ or both. The justification for these assumptions lies again 
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in the high-frequency behavior of the energy denominators of the Green's function G{k), 
analogous to the NRG assumptions for the self-energy. We can therefore use the same 
strategy as in the NRG approach for the self-energy: 

We again divide the full Matsubara frequency range into an L- and an if-region and, 
in the if-region, we calculate the i?-matrix at a higher temperature on a correspondingly 
coarser grid. For all required Matsubara frequency summations in the if-region, the R- 
matrix is then, again, interpolated onto the finer grid. In detail, this works as follows: 



The full i?-matrix, denoted by R^^^ at temperature Tj, is defined for matrix indices 



covering the full, dense icj^^-'-grid Ai, Eq. (^If), up to the maximum cutoff ^Iq. In all matrix 



multiplications during the i-th. NRG step, the full R^^^ is now replaced by a "diluted" R- 
matrix, which needs to be evaluated and stored only for indices on the NRG frequency 
grid LiU Hi. Those matrix elements which are eliminated by this procedure from the full 
i?-matrix are approximated by appropriate inter- and extrapolations from the i?*^*''-matrix 
elements retained. 

Consider, for example, a typical matrix-vector multiplication of the i?-matrix with a 
fermion frequency vector f{iuj). At the i-th NRG step, the required summation over all 
Matsubara frequency matrix indices icul^^ up to the cutoff, |ci;^*^| < is replaced by sum- 
mations over the "diluted" frequency grid L^U Hi. That is, rather than carrying out the 
matrix-vector multiplication on the full, dense Aj-grid to obtain 

g{^u;^^) ^ 5: i?»(.c.«,.a;«)/(^c.«), (40) 

n' 

we evaluate instead g{iuj'^^) on the NRG grid iu'^'^ G LiU Hi by linear interpolation of the 
summand in the ifj-region, which yields, similar to Eq. (pS]), 



i~l rp 
1 i' 
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= E ^:^^{^u^^.^u/^^)f{^iJ;!:^). (41) 

n' 

Note that the interpolation again introduces a temperature re- weighting factor Tji/Ti into 
the summation over the coarsened grid. The values of g{iuj^^) and /(iu;^) on the original full 
Aj-grid are then again representable by interpolation in terms of the g{iuj^'')- and f{iuj''^^)- 
values, respectively, on the NRG grid LiU Hi. 

This NRG grid representation of the i?-matrix, Eq. (|4T|) , is also used to carry out the 
matrix inversion of (5 = 1 + entering into Eqs. (|15D and (|16]). From Eq. (^I]) it is easy to 
see that, at the i-th NRG step, the problem is reduced to carrying out the matrix inversion 
of a "diluted" Q-matrix with matrix elements 

Q«(z4^-), ^o^ir^) = 5n^5,,, + , (42) 

where the matrix indices iu!^\iul^, ■* are restricted to the NRG grid LiU Hi. 

Based on the above-stated NRG assumption of approximate T-independence of R/T in 
the high-frequency region, we can express the i?^*^ matrix elements in the ^^Hi-Hi" region by 
i?-matrix elements already calculated in previous NRG steps. That is, analogous to Eq. (pHD, 
we have 

R^i^co^\^coi^P) = ^R^-"^ i^coi^\^coi^P) = ... = ^R^^"\^uJ\i\^uJ^?) (43) 

J-i-l -Lj" 

if both iuj^'^ and iool^, ■* G Hi. Here j" denotes the larger of j and j'. In other words, 
after the j-th NRG step, the Hj-Hj matrix elements of R^^^ are "frozen in" and need not 
be recalculated in subsequent NRG steps i > j, except for a change in the temperature 
prefactor. Only those matrix elements of R^'^^ connecting Lj to Hi and Lj to Li need to be 
calculated and updated to self-consistency during the self-consistency iteration of the i-th 
NRG step. 

In the pairing eigenvalue calculation, Eq. (0), we can use the same interpolative ap- 
proach described above for the i?-matrix, to carry out matrix multiplications with the pairing 
kernel, exploiting the diluted, but non-equidistant NRG frequency grid Lj U Hi at the final 
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temperature T = Ti. The non-equidistant nature of the grid will again introduce temperature 
re- weighting factors into the Matsubara summation, analogous to Eq. d^TI). 

The dimension of Q^^ on the NRG grid L^Uifi is Af{LiUHi)xAf{LiUHi) where Af{LiUHi) 
denotes the total number of grid points in Li U Hi. The inversion of Q^*^ at all momentum- 
energy-transfer vectors q [see Eqs. (P3|-P^)] is by far the most CPU time consumptive step at 
low temperatures. In, say, a "factor-2" NRG calculation, Af{Li U Hi) at low temperatures 
increases with 1/T as \og2{^^o/T); the amount of CPU time for inversion of a general 
DxD matrix scales as -D^; there are of order N^xNxNit such matrix inversions [one per 
g-vector and self-consistency iteration] in each NRG step, where is the spatial lattice (= 
k grid) size, A^^^ = Qq/{ttT) is the size of the original Matsubara frequency grid for maximal 
cutoff f^o and final temperature T, and A^it is a typical number of self-consistency iterations 
needed per NRG step; and the number of NRG steps to reach the final temperature T 
scales as \og2{-^^o/T). The total CPU time of the full NRG low-T self-energy calculation 
therefore scales as N^txNx (^Qq/T) x [log2(:^no/^)]^- This should be compared to the 
CPU time of a brute force calculation which, estimated along similar lines, scales at least 
as NitxNx (^Qq/T)'^. The savings in CPU time of the NRG approach relative to a brute 
force calculation is therefore a factor of order {^Qq/TY /[\og2{-^^o/T)]'^. 

The memory requirements for both NRG and brute force method are dominated by the 
storage of the i?^*^- and Q^*)-matrices. Following the foregoing dimensionality estimates, 
this scales as Nx{^no/T) x [log^i^^o/T)]'^ in the NRG and as A^x (i^o/T)^ in the brute 
force approach. Thus, in terms of memory consumption, the NRG saves a factor of order 
(^no/r)^/[log2(^fio/r)]^ relative to the brute force approach. 

IV. RESULTS 

To begin with, we explore the frequency and temperature dependence of the self-energy 
and Green's function of the pure Holstein model {U = 0) in Fig. ^, using a brute force 
calculation for temperatures T/t = 1.0, 0.5, 0.25, and 0.125. We use a maximum cutoff 
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frequency VLq/I = 25, a 16^ momentum grid, and an electron concentration (n) = 1.0. In 
Fig. ^(a) and (b), we plot the real part and imaginary part of the self-energy, respectively, 
as a function of Matsubara frequency at k = (1.77, 1.77), which is a point just above the 
Fermi surface along the diagonal direction. As expected in the high frequency region, the 
temperature dependence of the self-energy is much smaller than in the low frequency region. 
Note that we plot the imaginary part of as T,2/uJn instead of which emphasizes the 

low-iuj region at low temperature. 

Fig. ^(c) shows the Green's function at these parameters. It is quite clear that the 
high temperature behavior of the one-particle function is "frozen in" very quickly as the 
temperature decreases. As explained in the previous section, this is very important for the 
fermion frequency matrix NRG method to be applicable, since fluctuation propagators in 
the effective or pairing potentials, Eqs. ( p!^ - |T6| ) and (p3D, are calculated from the one-particle 
Green's function. Clearly, at high frequencies the Green's function can be interpolated by 
the high temperature results without loosing any informations. 

At Fig. ^ we plot the self-energy and Green's function right at the Fermi surface for 
k = (2.95, 0.20). At half-filling, the real parts of self-energy and of Green's function vanish 
at the Fermi surface, due to particle-hole symmetry. The low frequency behavior becomes 
much sharper at the Fermi surface, but the high frequency parts of S2 and G2 are almost 
temperature independent. 

In order to test the accuracy of our NRG approach, we use the factor-2 frequency NRG 
operation with constant Nj = 4 and Kj = 8, starting at To/t = 1.0, as described in the 
previous section. After three such NRG operation, we thus reach the final T = T3 = 0.125t. 
We plot the self-energy and Green's function for two different momentum points in Fig. |^, 
for k = (1.77, 1.77), and in Fig. ^ for k = (2.95,0.20). For comparison, we show the brute 
force results obtained for the same model parameters, temperature T, and Matsubara cutoff 
Qq = 25.13 (fixed). There is remarkable agreement between the two methods. However, 
the NRG grid Lj, Eq. (PB]), to be summed over in the self-consistency iterations, contains 
only 2Ni = 8 fermion Matsubara frequencies at each NRG step i, down to T = T3 = 0.125 
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and = 3.1. By contrast, in the brute force approach one has to sum over the full, dense 
Aj-grid, Eq. (p^, with Qq/IttT) = 64 fermion Matsubara frequencies. In table |I|, we list the 
memory and CPU time per iteration for both approaches. Savings of memory and CPU time 
of about two orders of magnitude are achieved by the NRG approach, without significant 
loss in accuracy. 

Next, we examine the stability of the NRG results against changes in the NRG control 
parameters and NRG protocol. Fig. |^ shows the s-wave eigenvalue calculation for three 
different "lower cutoff" frequencies Qe/t = 12.57, 6.283, and 3.142, which correspond to 
carrying out a total of, respectively, £ = 3, 4, and 5 factor-2 NRG steps, each starting from 
To/t = 4.0 and fio/t = 100.5, with fixed A^'q = A^i = ... = iV^ = 4. In each of these three 
calculations, the i-th factor-2 NRG step is followed by one "fixed-cutoff" step where the 
cutoff, fle+i, is left unchanged at fle+i = fit, while A^^+i is increased, beyond Ni, in order to 
lower the final temperature T = T^+i.lli This last, {£ + l)-th step is repeated with several 
different values of A^£+i, in order to scan the low-T regime. 

The s-wave is then determined by interpolation between two adjacent low-T points, 
T_|_ and T_, say, which bracket the instability, that is, their pairing eigenvalues, X{T^) < 1 < 
A(T_), bracket unity. The s-wave transition temperatures estimated from the three different 
calculations are Tc/t = 0.093, 0.090, and 0.091 for i = 3, 4, and 5, respectively. Thus, both 
the Tc and the A(T) results [Fig. |^ obtained with the three different lower cutoff protocols, 
£ = 3, 4, 5, are in excellent (~ 1 — 2%) agreement. 

For all further Tc results discussed below, the same "factor-2 plus fixed-fi^" NRG protocol 
is used, with an initial temperature Tg/t = 4.0, initial cutoff ^o/t = 100.5, and £ = 5 of 
factor-2 NRG steps, corresponding to Tf ft = 0.125 and a lower cutoff Qe/t = 3.142. 

Fig. p shows the s-wave as a function of the EP coupling strength Up. A 16^ grid is 
employed. The Einstein phonon frequency is Qp = l.Ot. Within the FLEX approximation, 
the s-wave Tc of the Holstein model increases with Up over a wide range of the EP coupling 
and saturates when the Up is comparable to the band width (8t). 

The finite size effect at the low temperatures needs to be treated carefully. We therefore 
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calculate the s-wave transition temperature as a function of electron concentration (n) in 
a wide filling range near half-filling. The results are shown in Fig. ^ for three different 
k-grids, = 16^, 32^, and 64^. The 16^ grid gives reasonable estimates over for a wide 
temperature range, down to the lowest temperatures we reach in the calculation, T ~ 0.02t. 
A 32^ grid, in general, covers this whole temperature range, down to 0.02t very accurately. 

Three sets of s-wave Tc are reported here corresponding to an on-site Coulomb U/t = 
(Holstein model), 2, and 3. All of them have the same Einstein phonon frequency Qp/t = 1.0 
and EP coupling Up/t = 4.0. The presence of the on-site Coulomb repulsion suppresses the 
s-wave pairing and Tc goes to zero or becomes smaller than our lowest numerically accessible 
T (~ 0.02t), when U approaches Up. Thus, in essence, goes to zero (or a numerically 
"very small" value) when the on-site Coulomb repulsion U overcomes the phonon-mediated 
on-site attraction Up. Note that there is no significant reduction of the s-wave-supressing U 
effect due to retardation, that is, due to the "pseudo-potential" reduction of the Coulomb 
repulsion. iiii This is perhaps not surprising, since, on the one-hand, the phonon-frequency is 
quite sizeable here, compared to the band- width 8t, and, on the other hand, near 1/2- filling 
there may be relevant electronic (charge, spin, and/or pair) fiucutation energy scales in the 
problem which are even closer to the phonon energy scale. Note also that, there exists an 
"optimal" doping, of about 20% to 30%, where a maximum s-wave transition temperature 
Tc occurs for this model in the FLEX approximation. 

A comparison between FLEX and conventional Eliashberg theory in the s-wave pairing 
regime U < Up, is shown in Fig. |10|. Here, we have used the bare potential U + Vp{iuJn — iuJn') 
as the effective exchange potential in the self-consistent (Migdal) self-energy calculation 
and as the pairing potential Vpair in the Ehashberg pairing eigenvalue calculations. Thus, 
in essence, our Eliashberg calculation negelects all the screening effects due to electronic 
particle-hole and particle-particle fluctuations which the FLEX approximation seeks to in- 
clude. At sufficiently large doping, > 15 — 20%, where converged results for Tc can be 
obtained in both approaches, the FLEX Tc is noticeably higher that the Eliashberg Tc. This 
suggests that the predominant {i.e., charge) fluctuations included in FLEX enhance the 
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s-wave pairing potential. Because of the lack of screening in the Eliashberg (and because 
of its presence in the FLEX) calculation, it is not surprising that the relative discrepancy 
between the two approaches becomes even larger in the presence of a finite on-site U = 2.0t, 
as shown by the dashed line results in Fig. 

A physically very interesting problem to study in the Holstein-Hubbard model is the 
competition between d- and s-wave singlet pairing. While the (1^2 _y2 pairing instability is 
by now well-established in several high-Tc cuprates, there is also mounting (although by no 
means unambiguous) experimental evidence for s-wave pairing in some of these materials.iii 
On the theoretical side, it is, within a weak-coupling self-consistent diagrammatic framework, 
well-established that a dx2_y2 pairing instability can be driven, or at least enhanced by AF SF 
exchange,iHi while being suppressed, due to self-energy effects, by the presence of phonon 
exchang e0a An s-wave instability, on the other hand, can be driven or enhanced by 
phonon exchange, while being suppressed by AF SF exchange and local Coulomb repulsions, 
as already discussed above. Thus the two possible candidate pairing mechanism for d- and 
s-wave pairing tend to be mutually destructive. 

The Holstein-Hubbard model, treated in the FLEX approximation, may be a reasonable 
starting point to investigate the magnitude these mutually destructive "anti-pairing" effects, 
both on the s-wave and on the d-wave side of the phase boundary. In Fig. ^(a) and (b), we 
plot the s- and, respectively, d-wave phase diagram at a fixed Up = 4.0t and, respectively, 
fixed U = 4.0t, for on-site Coulomb repulsions U/t = 0,2 and 3, and respectively, for on-site 
attraction Up/t = 0, 2 and 3. The Einstein phonon frequency is fixed at Qp/t = 0.5 for both 
cases. In comparing the two pairing states, we note that the s-wave state which has the 
higher Tc in the absence of its "anti-pairing" interactions {U = 0) is also suppressed more 
strongly, when its anti-pairing interaction U is turned on. We note also that the s-wave 
instability exists over a wider range of doping 1 — (n). An optimal doping with very broad 
Tc maxima is found in both s- and rf-wave pairing states. In the s-wave case this occurs at 
20% to 30% hole dopings, depending on the on-site Coulomb repulsion U/t; in the d-wave 
case the Tc-maximum is near 10% doping with a much less pronounced dependence on Up. 
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Another important feature to study in the Holstein-Hubbard model is the isotope effect, 
that is, the dependence of the superconducting on the isotopic mass M of the ions. In 
our parametrization of the model, this isotopic mass dependence enters only via the Einstein 
phonon frequency fip, that is, via Eq. (0), since all other parameters (t, f/, C, K, Up) are of 
purely electronic origin, i.e., not dependent on M. An important experimental measure of 
the lattice effects on Tc is the isotope exponent 

1 d\TiT, 



a = 



2d\nnp 



(44) 



ainM 

where the notation ...| means that the partial derivative should be taken with all above 



identified electronic parameters held constant. The second equality in follows from 
Eq. (D. 

In Fig. [1^, we plot both the s- and d-wave superconducting transition temperatures as a 
function of the phonon frequency Qp. In the s-wave case. Fig. |12|(a), Tc rises approximately 
linearly with Qp, up to Qp of about 2 — 3t. At larger Qp, Tc gradually becomes sub-linear 
and approaches saturation which is reached when Qp becomes of the order of the electronic 
bandwidth, that is, in physical terms, unrealistically large. The linear fip-dependence of the 
s-wave Tc implies that the isotope exponent is given essentially by its classical BCS value 
for conventional phonon-mediated s-wave superconductivity, 

a = i , (45) 

in the physically relevant low-Qp regime Qp 

In the d-wave case, Tc has only a very weak, slightly decreasing Qp dependence. The 
d-wave isotope exponent is therefore negative and and small in magnitude, typically with 

|a| < 0.05 (46) 

in the physically realistic f2p-regime Qp ^ t/4. This is much smaller than typical orders 
of magnitudes |a| ~ 0.4—1.0 observed in non-optimally doped cuprates, but confirms 
the conclusions from earlier calculations of the isotope exponent due to harmonic phonon 
exchange in diagrammatic d-wave pairing models.0'ii 
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V. SUMMARY 



In summary, we have developed an important generalization of the numerical renormal- 
ization group (NRG) technique for solving the self-consistent field equations of the fluc- 
tuation exchange (FLEX) approximation in the presence of a phonon-mediated, retarded 
bare interaction potential. In the presence of retarded bare interactions neither fast Fourier 
transform^ nor the previously developed NRG approachSIll for purely instantaneaous in- 
teractions can be employed. On the other hand, our generalized "fermion frequency matrix" 
NRG technique, produces large gains in computational efficiency, both in terms of CPU time 
and in terms of memory requirements, relative to a brute force calculational approach. 

In the physically most interesting low temperature regime, the CPU time and memory 
requirements of the brute force approach exceed by far the limits of currently available 
computational resources. By contrast, our generalized NRG method yield efficient, accurate 
solutions and allows detailed studies of superconducting instabilities in this regime, down to 
temperature scales 3 orders of magnitude below the electronic bandwidth. Our work also 
suggests possible avenues towards solving more complicated self-consistent approximation 
schemes, such the parquet approximationiS. 

We have tested and applied this NRG approach in the context of the FLEX approx- 
imation to the the 2D Holstein-Hubbard model. In this model, the FLEX equations are 
simplified due to the lack of momentum (q-) dependence in the bare Coulomb and in the 
bare phonon-mediated interaction potentials. However, more general electron-phonon inter- 
action models, including q-dependent potentials can easily be accommodated in our NRG ap- 
proach, by combining it with computationally efficient real-space basis representationsilBii. 
Such real-space basis representations have recently been implemented with great success to 
solve the FLEX equations in extended Hubbard model systems with instantaneous, but 
q-dependent bare interaction potentials. 

By varying the on-site Coulomb repulsion and EP coupling strength we have studied 
the competing s- and ci-wave superconducting instabilities of the 2D Holstein-Hubbard 
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model in the FLEX approximation. An optimal Tc for both cases was observed at about 
20% to 30% doping for s-wave pairing and 10% for ci-wave pairing. The s-wave phase 
is favored when the phonon-mediated on-site attraction Up exceeds the on-site Coulomb 
repulsion U and its transition temperature is suppressed by increasing U. Likewise, the d- 
wave phase is favored when the on-site Coulomb repulsion U exceeds the phonon-mediated 
on-site attraction Up and its transition temperature is suppressed by increasing Up. When 
U ~ Up, the Tc's of both pairing states are suppressed to zero or to a numerically inaccesible 
very-low temperature regime. 

Finally, the isotope exponent a for the s-wave state is BCS like, that is, a = ^, at 
realistic phonon frequencies Qp/t ^ 0.25. In the rf-wave state, the isotope exponents are 
negative and small in magnitude, with typically |a| < 0.05 in the physical phonon frequency 
regime Qp/t < 0.25. The overall magnitude of a is far too small to explain observed isotope 
data in non-optimally doped cuprates. Our full FLEX results thus support the conclusions 
of earlier d-wave isotope calculations by the present authors.0lllii 
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FIGURES 

FIG. 1. (a) The bare interactions in the Holstein-Hubbard model, which include the on-site 
Coulomb repulsion U and the EP coupling Vp{iiy,„). (b)-(d) show the bare vertices of the Hol- 
stein-Hubbard model in FLEX approximation. These diagrams include contributions (b) from 
density, and (c) from magnetic particle-hole fluctuations as well as (d) from singlet and 
from triplet particle-particle fluctuations. In (d), the upper (+) sign pertains to V^, the lower 
(— ) sign to y^. 

FIG. 2. Imaginary frequency discretization for the frequency NRG space, (a) Initial stage of 
the frequency space with cutoff J7o, including four positive and four negative fermion frequencies, 
corresponding to Nq = 4. (b) After one "factor-2" NRG operation, the lower cutoff is fli = 0,q/2. 
There are eight frequencies in "Li" and four frequencies in "ili". (c) After two "factor-2" NRG 
operations, the lower cutoff becomes = There are 8 frequencies in "L2" and eight 

frequencies in "^^2"- 

FIG. 3. Imaginary part of self-energy T,{k) = T,i{k) + iT,2{k) and Green's function G{k) for a 
brute force calculation at four different temperatures: T/t = 1.0 (cross symbols), 0.5 (squares), 
0.25 (circles), and 0.125 (lines) at k = (1.77,1.77). (a) Real part of self-energy Si. (b) Imaginary 
part of self-energy T,2/(jJn (c) The real part (solid symbols and line) and imaginary part (open 
symbols and dashed line) of G. The parameters are: = 25i, U/t = 0, and (n) = 1.0. 

FIG. 4. The same parameters as in Fig. 3 with the k = (2.95,0.20), which is right at the Fermi 
surface. Plots are shown for imaginary parts of self-energy and Green's function only; the real parts 
of self-energy and Green's function vanish on the Fermi surface at (n) = 1.0, due to particle-hole 
symmetry. 

FIG. 5. Comparison of the self-energy [(a) for real part and (b) for imaginary part] and Green's 
function [(c) solid line for real part and dashed lines for imaginary part] using factor-2 frequency 
NRG and a brute force approach. Results from 3 stages of NRG are represented by symbols. 
Results from brute force are represented by lines. The parameters are the same as Fig. 3. 
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FIG. 6. Comparison of (a) the self-energy and (b) Green's function using factor-2 frequency 
NRG and a brute force approach. Results from 3 stages of NRG are represented by symbols. 
Results from brute force are represented by lines. The parameters are the same as Fig. 4 

FIG. 7. Maximal eigenvalue As of the pairing kernel Eq. ( p2|) in the s-wave symmetry channel 
for three different lower cutoffs: O^/t =12.56 (solid line); 6.283 (cross symbols), and 3.142 (open 
circles). The model parameters are: VLp/t = 5.0, Up/t = 4.0, U /t = 2.0, and (n) = 0.75 

FIG. 8. The s-wave transition temperature Tc as a function of the EP coupling strength Up 
for the Holstein model {U = 0) with Einstein phonon frequency Vlp/t = 1.0 and electron filling 
factor (n) = 0.75. 

FIG. 9. Phase diagrams for the s-wave superconductivity of the Holstein-Hubbard model for 
different on-site Coulomb repulsion. The model parameters are Up/t = 4.0 and J7p/t = 1.0. 
The on-site Coulomb repulsion U /t = 0, 2.0, 3.0 (from the top curve to the bottom one). Three 
different k-grid sizes are employed. 

FIG. 10. Comparison of Eliashberg (open squares) and FLEX (open circles) solution for the 
phase diagrams for the s-wave superconductivity of Holstein-Hubbard model. The parameters are 
the same as Fig. 9 but data are shown only for = (solid lines) and 2 (dashed lines) . 

FIG. 11. Superconducting phase diagram of the Holstein-Hubbard model for (a) s-wave pairing 
with Up = 4.0t and several different U (b) d-wave with U = 4.0t and several different Up. 

FIG. 12. (a) s-wave and (b) d-wave transition temperatures Tc as functions of the Einstein 
phonon frequency Op near the optimal doping ((n) = 0.75 for s-wave and 0.90 for d-wave). In (a), 
Up/t = 4.0 and U /t = (solid line) and U/t = 2.0 (dashed line). In (b), U /t = 4.0 and Up/t = 2.0. 
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TABLES 



TABLE I. Comparison of memory and CPU time requirements between brute force and fre- 
quency NRG calculations. Parameters are the same as in Fig. 5. The calculations were performed 
on an IBM RS6000/397 workstation. 





memory 


CPU time per iteration 


brute force 


105 MB 


78.1 sec 


frofiuoncy NRG 


9 MB 


0.6 sec 
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